


Institutional Archive of the Naval Postgraduate School 


Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1980 


Analysis of combustion and heat transfer in a 
porous graphite medium. 


Vatikiotis, Costa Sozon 


Monterey, California. Naval Postgraduate School 


http://ndl.handle.net/10945/17515 


Downloaded from NPS Archive: Calhoun 


: Calhoun is the Naval Postgraduate School's public access digital repository for 
/ (8 D U DLEY research materials and institutional publications created by the NPS community. 
«ist : Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


NY KNOX appointed — and published -- scholarly author. 

; | LIBRARY Dudley Knox Library / Naval Postgraduate School 

411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 









i ee . ieee 










| , 8 a) ib é te 
| 4 Vi 
{ t my wrk aby te 
i] ea a! 
rh sito CP vena? 
‘ oa , 1, > a dh 
Wt Po a en Py th 
% rey Dal, ae ha " 
; . oF i % i js i t 
a0 1 any Daas t , 
ae “08 an ry hs | Ah fi 
sy BoE Cas Oe 
' Ses vv hic. oth 
ws ote Peved WT 
Soe th Bl i es 
vs 4 | Lee @ 
mf 2 a 
, “s } Ue . a) = 
. ’ . <a wer? A! * 
a *e ¢ % ‘ | 
. . ‘ : x rs 
_ 4 iy 
a 
_ | 
: = \ \ i] — i 
. ’ a’ 74 
\ : a Smo. 
a: Va 
im “4 ’ 
' ‘ ey 
‘ 7 ] t 
1e ev of pa 
4] , * wv * 
i + ' a 
A Ties. i 
» " ; ! 
aa | 
* 4 ¢ 
. 
1 hell | - 
* °° rn s Ait 
; . ar vl. /* i} 
ane 3 ; . 
° . . e '~ 
, = t = — - * 
RS i Ls — 
a = |= 
‘ - @e 
4 ° » 
° - -_ — = + = 
7. = 
° ‘ ~ - ‘all i} -—_ 
4 - zi - 
’ a ' : 
“ * a 
. 3 = 
am Ss =, 
. —— 
a z - -— 
™» . -_ = 
= a eg sen — 
~ ; . -— 
. © L . 
- - 4, — - q 
F ii = —— 
<- 4. s 
4 - - = “= 
w = 
aay ' 4 
> 
a> oe 
. P ‘ - - 
. = . 
mit 
F Ny 
ry 
e ii 
! . 
i) 
' - ~ 
” Pe 7 re 
> 
* sd i, a 
' di @ 
Us Ty ar q 
| = so 
= 
re 
° é 
S ’ 
¢ 
é 
e 
° 3 a L] 
é 
re] : | 2 
i] : . 
r ¢ é ; 
° 7] ' 
i] 
. : é . 
. / ' 6 
? it ° ) * 
3 
6 
' ‘ , 
i e 
’ f ' 
. } . r i} : iy 
- ‘ , b ‘ 
% i i Th ’ ay | Lt fa al 
! ae C0) MS eee 
Mi P , ac? x 
; | 7 | ; ' | = a, ; a 
, fe th , i ‘i ' iti tee, LF ins, J 
r , ‘ 5 \ \ ‘i Dy eS 1) i. CARA 


DUDLEY KNOX LIBRARY 
NAVA. POSTGRADUATE sCHOOL 


MONTEREY CALIF 93946 




















NAVAL POSTGRADUATE SCHOOL 


Monterey, California 





BHESIS 


ANALYSIS OF COMBUSTION AND HEAT TRANSFER 
IN A POROUS GRAPHITE MEDIUM 


by 
Costa Sozon Vatikiotis 


June 1980 


Thesas AGVIUSOr: D. Salinas 





Approved for public release; distribution unlimited. 


7196285 








UNCLASSIFIED 


SECURITY CLASSIFICATION OF THIS PAGE (Wren Dare Entered) 


REPORT DOCUMENTATION PAGE 


2. GOVT ACCESSION NO 





READ INSTRUCTIONS 
BEFORE COMPLETING FORM 


5. TYPE OF REPORT & PERMIOO COVERED 
Engineer's Thesis; 


June 1980 


6. PERFORMING ORG. REPORT NUMBER 


















TITLE (and Subsettte) 
Analysis of Combustion and Heat Transfer 


in a Porous Graphite Medium 















. AUTHOR 4) 8. CONTRACT OR GRANT NUMBER(e) 






Costa Sozon Vatikiotis 











9. PERFORMING ORGANIZATION NAME ANDO AOORESS 
Naval Postgraduate School 


Monterey, California 93940 









CONTROLLING OF FIZTE NAME AND ADDRESS 


Naval Postgraduate School 
Monterey, California 93940 





12. REPORT CATE 
June 1980 
MONITORING AGENCY NAME & BODOMESE( II ailferant (rem Contreiling Office) 1S. SECURITY CLASS. (of tate report) W 


Unclassified 


Sea. DECL ASSIFICATION/ COWNGRAOING 
SCHEOULE 


Approved for public release; distribution unlimited. 






OISTRIBUTION STATEMENT (of thie Repert) 





DISTRIBUTION BTATEMENT (fal the aharcect entered in Bieck 20, ({ dilferent trem Report) 


. SUPPLEMENTARY NOTES 





. KEY WORDS (Continue on reverse side ii neceseary and identify vy block number) 
Combustion, Heat Transfer, Carbon, Graphite, Fiber, 
Porous Media, Conduction, Convection, Radiation, Diffusion, 
Exothermic Reaction, Nonvolatile, Finite Element Method, 
Galerkin Formulation 









ABSTRACT (Continue an reveree side if neceeeary and identify by bleck number) 


The problem of a porous graphite fiber plate subject to 
combustion is formulated. The transient one-dimensional model 
leads to two heat transfer equations on the graphite and air, and 
amass transfer equation on the oxygen. In addition to a combustion 
heat generation term, the heat transfer model includes the 
mechanisms of conduction, radiation between fibers, and convection 
due to induced air flow through the mat. The mass transfer model 


20. 









OD . en 73 1473. eaition oF | wov 681s omsovere 


(Page 1) S/N 0102-014- 6601 | 
] SECURITY CLASSIFICATION OF THIS PAGE (Prem Dete Entered) 


UNCLASSIFIED 





GCOcumMry CLASSIFICATION OF TwIS PAGE(When More Entered. 


DD Form. 1473 
R ‘ah ( 
S/ 0102-014-6601 2 STCUAITY CLABSIFIC ATION QF TiS PAGESWRen Date Eniored) 





#20 - ABSTRACT - CONTINUED 


includes diffusion and convection mechanisms, as well as a 
combustion consumption term. The temperature dependency of 
the system parameters is taken into account. The resulting 
nonlinear, coupled partial differential equations are solved 
by a Galerkin formulation of the finite element method. A 
number of computer analysis results are presented. The 
results show the effects of initial conditions, plate thick- 
ness, fiber diameter and air flow rate on ignition and 
extinction. 


UNCLASSIFIED 





Approved for public release; distribution unlimited. 


Analysis of Combustion and Heat Transfer 
in a Porous Graphite Medium 


by 
Costa Sozon Wark TOrrs 


Lieutenant, United States Navy 
bPaompeeleoriaa@=atate University, 1971 


Submitted in partial fulfillment of the 
requirements for the degrees of 
MASTER OF SCIENCE IN MECHANICAL ENGINEERING 
and 


MECHANICAL ENGINEER 


from the 


NAVAL POSTGRADUATE SCHOOL 
June 1980 





DUDLEY KNOX LIBHA! 
NAVAL POSTGRADUATE SCHOU 


MONTEREY, CALIF 93840 
ABSTRACT 


The problem of a porous graphite fiber plate subject to 
combustion is formulated. The transient one-dimensional 
model leads to two heat transfer equations on the graphite 
and air, and a mass transfer equation on the oxygen. In 
addition to a combustion heat generation term, the heat 
transfer model includes the mechanisms of conduction, radia- 
tion between fibers, and convection due to induced air flow 
through the mat. The mass transfer model includes diffusion 
and convection mechanisms, as well as a combustion Son gumiseden 
term. The temperature dependency of the system parameters 
1s taken into account. The resulting nonlinear, coupled 
partial differential equations are solved by a Galerkin formu- 
lation of the finite element method. A number of computer 
analysis results are presented. The results show the effects 
of initial conditions, plate thickness, fiber diameter and 


air flow rate on ignition and extinction. 
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im INTRODUCTION 


The combustion of a composite laminate plate consisting 
of graphite fibers in an epoxy matrix is discussed. The 
combustion process for a graphite laminate takes place in 
two stages, epoxy combustion followed by combustion of the 
graphite. At the end of the first stage the epoxy has burned 
away exposing a porous graphite mat. Spacing between the 
fibers 1s maintained by the residue of the combustion pro- 
ducts. This present work is concerned with the last stage 
of the combustion process, specifically, the thermal response 
of the remaining graphite fiber mat. 

The selection of the second stage for the analysis was 
made as a result of observations during composite plate 
"burn" experiments at the Naval Weapons Center, China Lake [1]. 
The objective of the experiments was to assess damage of com- 
posite aircraft structures subjected to open deck fires on 
aircraft carriers. For the experiments, an epoxy-graphite 
laminate was mounted into the wall of a wind tunnel parallel 
to the flow. One side of the composite plate was exposed to 
the wind tunnel flow and the opposite side was open to the 
environment. This arrangement was to simulate a composite 
wing section or fuselage subjected to typical wind condi- 
tions present on a carrier flight deck. It was observed 


that after being subjected to a fire which resulted in the 
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epoxy burning away, the graphite fibers would continue to 
smolder or burn. In other instances, the fibers would cool 
depending on the flow velocity and plate thickness. To 
understand this behavior better, a mathematical model was 
formulated to predict the conditions for which the exo- 
thermic process is self-sustaining. 

The air flow over one surface produces a pressure 
differential across the plate which induces a convective 
air flow through the porous medium governed by Darcy's Law. 
As a result, there is an enhancement of internal convection 
heat transfer, as well as a source of oxygen for combustion. 
The heat transfer mechanisms included in the model are 
(1) conduction, (2) convection, and (3) radiation. In addi- 
tion, non-volatile combustion is included in the energy 
balance as a heat generation term of Arrhenius type. The 
oxygen concentration is accounted for by a molecule-mass 
balance which includes (1) molecular diffusion, and (2) mass 
transfer by convection. Consumption of oxygen due to combus- 
tion is accounted for by a term similar to the heat generation 
term. 

The formulation of a one-dimensional model is presented. 
All properties are treated as temperature dependent in the 
transient analysis. To account for heat transfer mechanisms 
between the air and porous graphite medium, the air and 
graphite temperatures are treated as independent variables. 


The energy balance on the graphite and on the air, and the 
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molecule-mass balance on oxygen result in a system of three 

nonlinear, coupled partial differential equations. 
Integration of the field equations is accomplished by 

a Galerkin formulation of the finite element method. Due 

to the inherent stiffness of the field equations in the time 

domain, a modified implicit-Gear integration scheme developed 


by Franke [2] is used for a more efficient solution. 
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II. THEORY AND BACKGROUND 


A. INTERACTION OF HEAT TRANSFER AND COMBUSTION 

In this work, the combustion model proposed by Semenov 
was adopted. It is described in the texts of Frank-Kamenetskil 
[3] and Vulis [4]. A brief discussion of those features of 
the model which relate to the present investigation will be 
given. Fundamental to the model is the relation of reaction 
rate to temperature, and the interaction of the heat genera- 
tion and the heat transfer from the system. In general, the 


reaction rate, r, may be shown by the Arrhenius law as, 
r = A(T,¢) exp(-E/RT) 


where ies is the characteristic time of the chemical reac- 
tion, E is the activation energy, R is the universal gas 
constant, T 1s absolute temperature and ¢ 1s the concentra- 
tion of oxygen. The concentration appears as oe nel A(T, 4) 
for an nth order reaction. The heat generated by the exo- 
thermic process is obtained by multiplying the reaction 
rate, r, by the enthalpy of formation of the combustion 
Products. 

In Figure 1, the heat generation, eee 1s plotted as a 
function of temperature. The curve is referred to as the 


S-curve for apparent reasons. The S-curves are distinguished 


1S, 





by two regions. Lower temperatures and reaction rates 
characterize region I, while region II is characterized by 
higher temperatures and reaction rates. In region I, for 
low reaction rates, there is an abundant supply of oxygen. 
The reaction is controlled by the temperature and the region 
is known as the kinetic regime. Here, the reaction rate 
increases exponentially with increasing temperature. In 
region II, the reaction is limited by oxygen concentration 
and is known as the diffusion regime. The reaction's weak 
dependence on temperature in region II is shown by the char- 
mereristrc flattening Of the S-curve in Figure 1. For higher 
concentrations of oxygen, the diffusion region is associated 
with higher reaction rates. The interaction of heat genera- 
ero , Sa and heat loss, Gos will now be described. 

In addition to the S-curve, Figure 1 shows two heat 
loss lines, Gol and Goo° They represent the heat transfer 


by convection for air temperatures, T, and T respectively. 


2 


The intersection of the heat generation curve and the heat 


loss curve, oe =? defines the stable gquasi-stationary point, 


A. As the air flow temperature increases from T,; to T the 


Qe 
quasi-stationary temperature increases continuously from Ta 


to Tr: FOE cemperacure Th! Goo 


infinitesimal increase in the air flow temperature results 


is tangent to ac and an 


in a jump of the reaction temperature from T_ to T At point 


I Be 
ieevwaten 2s detinea as the “critical ignition condition", 


the reaction moves from the kinetic to the diffusion regime. 
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Temperature T. is referred to as the ignition temperature. 


I 
Frank~Kamenetskii states the ignition temperature of graphite 
iS approximately 1100 to 1300 degrees Celsius. In order to 
keep the reaction model simple, the present investigation was 
limited to the kinetic regime. Although the description of 
the heat transfer process just presented is for convection 
only, the underlying ideas are valid for the conduction and 
radiation heat transfer mechanisms as well. The slope of the 
heat loss lines, qos is equal to the product of the internal 
heat transfer coefficient and the surface area of fiber per 
unit volume. The slopes would not appear constant, and slight 
curvature in the heat loss lines, Gos is apparent when radia-~ 
tion effects are included. As shown in Figure l, an increase 
in the heat transfer rate resulting, for example, by an 
increase in the internal flow rate would result in a lower 
graphite temperature for a given air temperature. 

One reason for limiting the model to the kinetic regime 
was due to the jump in temperature and reaction rate at 
ignition. Though the discontinuity may be accounted for, the 
results show that for the region of temperatures approaching 
the diffusion regime, the heat generation will dominate the 
process. The system temperatures then rise rapidly. Other 
factors that limit the model to the lower temperatures are 
(1) the chemical reaction fixing the heat of combustion (com- 
plex at higher temperatures) and (2) the rapid consumption of 


fibers after ignition. 
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Eee HEAT TRANSBERSEQUATIONS FOR POROUS MEDIA 

In their investigation of oil recovery from underground 
reservoirs, Green and Perry [5] developed heat transfer 
equations for the solid and fluid phases in a porous medium. 
Their model included the three basic mechanisms: (1) physical 
movement of the fluid which carries its own heat capacity, 
(2) conduction of heat through the solid and fluid phases, 
and (3) convective heat transfer between the solid and fluid 
phases. Radiation heat transfer was assumed to be negligible, 
and combustion was not a consideration. The differential 


equations for the fluid and solid phases 


2 








oT; aT e) Te 
eeu C : PeCeP aE = TU eC EP Ot K ep 2 = hz(?,. - TQ) lk L) 
oT. son, 
Solid: pC, (l-PIZE = ki (1-p) a: oe hz(T, - T.) CEEs 2) 


were solved numerically by the finite difference method. 

In equations II.1 and II.2, c is the specific heat at constant 
pressure, u is the fluid velocity through the medium, h is 

the coefficient of heat transfer, and k 1s a pseudothermal 
conductivity. The model did not include change of properties 
with temperature. Their numerical results showed good agree- 
ment with their experimental results. Two useful conclusions 
derived from their analysis were: (1) both conduction and 
convection heat transfer mechanisms are important for internal 


Reynolds numbers less than one, and (2) for values of the 
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dimensionless parameter, & = (h2/k,)°/* (k,/0 pc eu) greater 
than .342, the fluid temperature is approximately the solid 
temperature. Riaz [6] proposed the following equation 


resulting from the equal fluid-solid temperature assumption, 


aT aT 
Psspe * PeCe4Gx = *g 7 =) 
where the coefficients are defined as those in equations 
II.l and II.2. For the present model, the more general 
approach (1.e., Green and Perry) will be taken in the formula- 


tion of the energy equations. 


aE 





The coORMUEATTON OF THE FIELD EQUATIONS 


A. INTERNAL FLOW MODEL 

Referring to Figure 2, the porous mat may be represented 
by a flat plate. One side is exposed to still air, and the 
Seeger tO an air flow or velocity, U_. The environment on 
both sides is at ambient temperature, ess 

The first step in developing the model was to derive 
an expression for the flow through the porous plate. Effects 
of the combustion process on the physical properties were 
neglected. Specifically, these were (1) variations in the 
porosity due to fiber consumption, and (2) effects on the 
viscosity and density of air from the introduction of combus- 
tion by-products into the flow. Preliminary calculations 
showed that the porosity and density of the material did not 
change appreciably until temperatures reached approximately 
2000 degrees Fahrenheit. Further, since eighty percent of 
air by volume is essentially inert, it was reasonable that 
its properties would not be significantly altered during the 
combustion process. Assuming a heterogeneous material of 
constant porosity, a Reynold's number for the interior flow 


is defined as 


Re. = pu ,d/u (oeiiies 1) 
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where O. is the mass density of air, Be is the pore velocity, 
dis the diameter of the fiber, and u is the dynamic viscosity. 
Extensive experimental work shows that Darcy's law for flow 
through a porous medium is valid for a limited range of 
Reynolds numbers. Muskat [7] proposed that Darcy's flow is 
valid for Re, < 1. Scheidegger [8] points out that a number 
of investigators give a much higher upper limit. Preliminary 
calculations of Re. for the model showed typical magnitudes 

On the order of 1. Neglecting the gravity term, Darcy's law 


for flow in porous media is, 


_ m dP 
4, = 7 GE CETL 2 ) 
assuming a constant pressure gradient, equation III.2 
becomes, 
ligl (SLE 
= - =- = Sg ees 
U, a ( ) 


where AP is the pressure differential across the plate 
thickness, m is the specific permeability of the mediun, 
and u is the dynamic viscosity. In equations III.1 and III.3, 
Ps and pup were treated as temperature dependent. 

Specific permeability is a measure of a porous medium's 
ferlicy co allow £luid to flow through. It is dependent only 


on geometry and the physical nature of the medium. It has 


dimensions of length squared. Specific permeability is often 


Pa Mt 





called the hydraulic conductivity due to the similarity of 
Darcy's law to that of Fourier's law for heat conduction. 

The permeability for a particular medium is normally measured 
by experimentation. However, there are several empirical 
models that may be used to obtain values for permeability that 
are in agreement with experimental results. For the idealized 
geometry of the porous medium shown in Figure 2, the permea- 
bility was obtained by a serial type model, equation III.4, 


proposed by Scheidegger [8], 
Il 2 


where the porosity, p, is defined for a prous material 


comprised of fibers with cylindrical shape as, 


Diaz lee 7 (a/D) * (III.5) 
and D is the thickness of a ply. Porosity has units of 
volume of void per unit volume of medium. The average pore 
diameter, § was defined as, 

6 = (s + D)/2 
for the geometry in Figure 2. The tortuosity, t, 1S a measure 


of length of travel for a fluid particle per unit thickness 


of the medium. Referring to the geometry of Figure 2, the 


De 





tortuosity depends on the ratio, dad/D. The lower limit 
occurs when d/D equals zero and for the upper limit, d/D 
equals 1. Carman [9] presents a table of tortuosities for 
given geometries and particle shapes. For the model, 
Carman's suggested value of t = 1.4 was used. 

To obtain the pressure differential, AP, across the 
porous plate, Bernoulli's equation was used. The following 
observations showed this to be a valid assumption. Schlichting 
[10] points out that for the ratio of nay Os in the range of 
meoor to .01, the effects of "blowing" or "suction" on the 
potential flow over the porous plate may be neglected. A 
Eypical value of ae for the model at which U_ = 25 kt. was 
.0028. For steady flow over a flat plate, the flow field 
outside the boundary layer may be described by Bernoulli's 
equation. This is a direct result of the Navier-Stokes equa- 
tion. In addition, the pressure gradient across the boundary 
layer may be taken equal to zero. Therefore, Bernoulli's 
equation (eq. III.6) may be used to obtain the pressure 


differential across the plate, 


2 
f= ax + = constant. (III.6) 


2 
a 
The parameters in equation III.6 are defined as follows: 
P, pressure; as the density of air; U, velocity of the air 
over the flat plate. For the model, the density of air was 


approximated by the Ideal Gas law as 


Zo 


rh - 





as P/RT. Camis 


“Aw 


where Re 1s the gas constant for air, T. is the absolute 
temperature of air, and P is the pressure. Substituting 


this into equation III.6, gives 


nN 


Ree 2 


inf ===) ap + = = constant 





and upon integrating, the expression becomes 


5 2 
gg BN Jog) Ee CS constant 
aa 2 
Si, 
v2 v2 
Letting 
eee Ut, = 0, T., 21,5 = Ter Pp = Ppe Un = Uy 
and substituting these in above, yields 
Oia 
= ie 
RT, inP, RT inP, 5 


Rearranging the previous expression gives, 


24 





| 2 
son ee IB! (P/P,) = =F 


Taking the exponential of both sides and solving for PL 


mesults in 
2 aA 
Bo = See (Sy aa 


L 


From the above expression and noting that AP = P, -P, 


AP may be expressed as 
2 “a 
AP = Plexp(-U,/2R_T,) - 1] 


Darcy's law in approximate form (eq. III.3) for the model 


is 
and substituting for AP, it follows that, 


Pp 
= tl | oe oe” 2 
u = i stl exp ( U/2RT,)]} (SEICIE 6183), 


where 
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The pore velocity, aoe of equation III.8 is for an 
isothermal medium at ambient temperature, Dy . Since the 
model is to investigate the transient problem, a general 
expression to obtain pore velocity at different temperatures 
had to be developed. 

Consideration of the continuity equation for one-dimen- 


sional steady flow, 


d(o u,,) /dx a) (ied!) 
becomes 
ur tae 7 ox) + p , (du /dx) = Q 
‘ene 
EL, Me ayiees a De Mere CES) 


Multiplying through both sides by dx, the expression 


becomes, 


ud = du 
a Weve ae 


Separating variables and integrating both sides, gives 
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Taking the exponential of both sides, gives 


Rearranging, the pore velocity is expressed as 


u _ u GAL ese 6) 
p P., yer! Pa ( 
Referring to the Ideal Gas law for the density of air, and 


substituting this into equation III.10 gives 


Le = (0 Uy eRaT,)/P Cite. 1) 
The pore velocity is now a function of air temperature and 
pressure. However, in preliminary calculations, the pressure 
difference, AP was small and for the range of parameters of 
interest, it did not exceed five percent of the magnitude of 
ambient pressure. Therefore, as a simplification, it was 


assumed that P equals the ambient pressure in equation III.1ll, 


ae 





and the expression becomes 


aA 


ae “= aloe aa 
Noting that op. = ey eee er the pore velocity may be expressed 


as 


u Seed T/T (Tae. 1 2 ) 
where bee is obtained from equation III.8. a angie, aise 
the absolute temperatures of the air at a point in the interior 


of the medium and of the environment, respectively. 


B. ENERGY EQUATIONS 

Here we consider the development of the energy equations. 
However, before continuing, discussion must be made as to 
the approach taken for energy balances in porous media. [In 
previous work by Denbigh and Turner [ll], and by Colladay 
and Stepka [12], it was suggested that the temperature of 
the porous solid and of the fluid be considered equal. This 
greatly simplifies the formulation of the problem. However, 
under certain conditions, as shown by Green and Perry [5], 
the assumption of equal temperatures may not be valid. With 
recent developments in numerical techniques permitting the 
solution of non-linear systems, the development of the present 


model for heat transfer in a porous medium was not restricted 
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to the assumption of equal temperatures. Therefore, the 
temperatures of the air and of the porous solid are treated 
as independent variables in the present formulation. 

The model is based on two assumptions. First, due to 
the small thickness of the graphite fibers, the temperature 
across each individual fiber was assumed to be constant. 
Second, for a similar size consideration, the temperature of 
the air in each individual pore was assumed constant. How- 
ever, the temperature from fiber to fiber and from pore to 
pore was not restricted, and was allowed to vary according 
to the heat transfer mechanisms described next. 

To perform energy balances on both the porous solid and 
on the air, a differential volume of porous material was 
segregated into respective volumes of constituents, that is, 
graphite fibers and air (Figure 3). An energy balance may 
be done on each volume separately. The convention used for 


the balance is: 


Increase in 


Heat into Heat _ Heat out 
av * Generation ~~ of av Internal 
Energy 
where 
Ovo = dx dy “dz 
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ly Fiber Heat Transfer Equation 


Considering the differential volume, (l-p)dvV of 
graphite fibers, and "smearing" the fibers to form a macro- 
scopically homogeneous material, the heat transfer mechanisms 
are shown for the x-direction in Figure 4. The heat transfer 


mechanisms are: q heat conduction through the medium; 


cond’ 


q , convection from the fibers to the air; dead’ radiation 


Conv 


transfer from fiber to fiber; dur (To1¢) heat generation rate 
per unit volume. AH is the enthalpy of formation, and 
a7!) is the fiber mass consumption rate. This term will 
be discussed in detail later. 

Representing the terms in Figure 4 by Taylor series 
expansions and combining them with respect to the convention 


previously stated, the energy balance is 


ere | 
cond 
+ aur (T,,9) (1-p) dxdydz doond ax 


doond a Trad 


9q 
+ q 4 = ex + ¢ 


JE 
rad * “x eee 


conv 


Subtracting terms, q and dea from each side and 


d 


rearranging, the expression becomes (neglecting the higher 


cond 


order terms) 





ret aq 
cond rad . 
— x or al a8 OO) Me p) dxdydz 
= 22 (1-p) dxdyd Gea 1s 
= op (l-p) dxdydz Ssieo) 


30 





The heat transfer terms in equation III.13 will now be 
considered. 
a. Thermal Conduction 
Fourier's law for conduction may be written for 


the incremental volume as 


= Kya (III.14) 
where He is the effective conductivity of the porous solid, 
and ‘e is the temperature of the fibers. There are a variety 
of models for the effective conductivities of porous media. 
They depend largely on the geometry and on the nature of 
the constituents. For the cylindrical shape of the fibers, 


the empirical relation proposed by Rohsenow [13] was used 


2 72 
k = ke {ee tl ns a ee my ) 


eet 5 15) 
J fea. M 


for A < 1. The parameters of equation III.15 are defined 
as follows: 8 = d/D (refer to Figure 2), y is 

1 A 

y = Alay] 

2 (1l=A) 
where A = be kK. is the conductivity of the air; and Re 
is the actual conductivity of the graphite in bulk form. 
Since Kk. depends on temperature, Ne was also treated as 


temperature dependent. 


ek 





b. Radiation Heat Transfer in Fibers 
Radiation heat transfer in the model was repre- 


sented by an analog to Fourier's law of conduction as 


oe 


ae _ 9 
ame k. = ay dz Grant. 16) 


where Kk. is an equivalent conductivity of the fibers due 

to radiation. Simple assumptions and approximations may be 
made to show this. Assuming air to be transparent to radi- 
ation, and treating a ply of graphite fibers as an infinite 


wall, the net heat flux between plies may be written as 


_ Oe (Tr? = 74 


rad (Zee) ax gxtdx? (ieee 7 ) 


q 


where Tox and Tox+dx are the respective absolute ply tempera- 


tures, « is the emissivity of the graphite, and o is the 
stefan-Boltzman constant. Treating ee as the variable 


temperature, q may be expanded ina Taylor series about 


e 
rad 


eredx" Neglecting the higher order terms, the series 


expansion may be written as 


oe ,24 aa age 3 if a 


i ( ee eed ae Toxtax!T gx 7 Tgxtdx’ 


dead 2-€ 
Simplifying, the above equation becomes, 


- —_ 4oe 23 a me 
q 7 Nope cbe eaes T oxtdx 





) (ite 1es)) 


a2 





Taking Fourier's law for heat transfer, 


aT 
dead =) aaa Se ax 
oul 
and approximating = with 
qT _ ong _ T oxtdx 7 T.. 
ao.” ax AX 


See becomes 


”~ aw 


ah ee gxtax 7 Fgx! 


rad Ve AX (ITT.19) 


Multiplying the right side of equation III.18 by Ax/Ax, 
equation III.18 and equation III.19 may be equated 


“aw “aw wN w~N 


= St ek. 
" a, 5 Tgxtdx T gx pe esac Ax 73 gxtdx gx’ 
drad r AX Zoe gx+dx AX 


Cancelling terms, Kk. becomes 


Seek x. ,3 
ie 2-€ gxtax 





lees 2:0} 


where AX iS now equal to 6, the average pore diameter. 
From the close spacing of the fiber layers, the temperature 
difference will be small as compared to the magnitude of the 


temperature. Noting this the average absolute temperature of 


a) 





the fibers, . may be substituted for T The equivalent 


gGxtdx~ 


radiation conductivity expression becomes 
Crerr. 21) 


and the radiation heat transfer from fiber to fiber may be 


represented by, 





—3 dy dz tee 22) 


If the temperature gradient in the fibers is not 
large, then equation III.22 will be a good approximation of 
the radiation heat transfer between the fibers. 

c. Convection Heat Transfer 

The convection heat transfer term in Figure 4 
was treated in a similar manner to that of the one-dimensional 
fin equation. gq was introduced as 


CONV 


ee = hy(T,-T,) dA (iit.23) 
where Be 1s the temperature of the graphite, To is the 
temperature of the air, dA is the surface area of graphite 

in the differential volume, and h; is the internal convection 
heat transfer coefficient. An empirical expression in the 
form of a Colburn j-factor was developed by Yoshida, 


Ramaswami, Hougen [14]. This was used to determine he, 
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Dy a") 2/3 
c,G k 
bo a 


= 9] Bee ai (III.24) 





mom Re < 50. Cy 


medium and is defined as, 


is the effective specific heat of the porous 


oe PS. + (Topic. 

where Cc. PmioesveerrTeCeneat Of alr at constant pressure; 
os is the specific heat of graphite; and p is the porosity 
as defined by equation III.5. Ga is a pseudo mass velocity 
defined as Go = Pate. Also in equation III.24, uw is the 
viscosity of air, and kK PS eENemconauctivity Of alr. The 


Reynolds number appearing in equation III.24 is defined as 
Re = G/2 u a 


and is not the same as Re; defined previously for Darcy's 

law. zis the surface area of graphite fibers per unit 

volume, and is determined from geometrical considerations 

which will be described shortly. The parameter, Yo? is a 
dimensionless shape factor which depends on the geometry of 

the fibers. For a cylindrical fiber shape, ie is equal to 

.91. From the assumption of constant porosity, z remains 
constant. However, the air properties are temperature dependent; 


therefore, h; is allowed to vary with temperature. 
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d. Heat Generation Rate 
The reaction rate, BA whether based upon 
the Collision Theory or on the Theory of Absolute Reaction 


Rates, results in an expression of Arrhenius type, 
r = A(T,¢) exp(-E/RT) 


where E is the activation energy, R is the universal gas 
constant, and T is absolute temperature. Although E is well 
defined for the combustion of graphite, A(T, $) is a function 
of temperature and concentration and depends upon chemical 
kinetic theory. Development of a general model for the chemi- 
cal kinetics is beyond the scope of this work. A less ele- 
gant, but useful approach of using a relation obtained 
experimentally was taken. A relation derived from the work 
of Parker and Hottel [15] was used to predict the combustion 


rate of graphite fibers, 


- Apcitegao’ R. te’? y 6 exp(—ees83) —R (rir.25) 
g T, ft~-hr 


The development of equation III.25 from the original work 
may be found in the appendix. This expression assumes a 
Simple first-order reaction for the combustion of graphite 
in air. However, Frank-Kamenetskii [3] has further refined 
Parker and Hottel's expression to yield for the present 


model, 


6 





eo x 10  y 6” exp ( 


pS) ic igeys) lbm 
na ee a 


) 
ab ee ae 
g 


(iio: 


The combustion is now characterized as an nth order react- 
tion, where n is in the range of 1/3 to 2/3. This form of 
the reaction rate expression better approximates Parker and 
Hottel's experimental data. For the analysis, n was chosen 
to be 1/2. However, it should be pointed out that the behavior 
of the reaction varies significantly within the range of n 
proposed by Frank-Kamenetskill. 

In limiting the model to the kinetic regime of 
combustion, it was assumed that the dominant chemical reac- 


mon is 


Ca 0 CO 


This is an idealized treatment of the reaction that actually 
takes place. Kanury [16] states that, the actual reaction 
1S more complex, yielding various concentrations of carbon 
dioxide and carbon monoxide, depending on the temperature. 
Noting the abundant oxygen supply present in the kinetic 
regime, it is reasonable that the leaner oxide is the preva- 
lent by-product of combustion. The fuel-oxygen ratio, f, 
for this mechanism is 12/32. Carbon monoxide may be formed 
when there is little free-oxygen present, such as, during 


combustion in the diffusion regime. The stoichiometric 
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fuel-oxygen ratio, f = 12/16, for this mechanism is naturally 
greater. To obtain the heat generation rate, the reaction 
rate, r, was multiplied by the enthalpy of formation, AH, 
of carbon dioxide. AH has units of Btu per lbm of graphite 
consumed. 

e. Change in Internal Energy 


The last term of the graphite energy equation, 


the rate of change of internal energy, <=, results by 
setting E = ce T. hus 
g Pg gg , 
.E oan 
xp (L-p) dxdydz = sO) Pe Te dxdydz Ciniedien 27 ) 


where 2s 1s the density of graphite in bulk; Se is the 
specific heat of graphite; and p is the porosity. 

The heat transfer mechanisms for the graphite 
Peeers, equations IIIT.14, I1It.16, III.23, III.26, IIrI.27 are 
substituted into equation III.13 which upon dividing through 


by dV = dxdydz yields the heat transfer equation for the 


fibers, 
aT oT 
a _ go ap) oe my GA pm _ 
(1-p) ax (Kg al + 1 Eine sa hs ag ae co 
oT 
+ ea = (Lope oc, — 


The coefficient, (l-p), accounts for the fact that not all 


of the volume of porous medium is comprised of graphite. 


ae 








Pacing through by (l-p), combining the first two terms, 
and defining dA/dV as z, the surface area of graphite per 


unit volume, the expression becomes 


: oT hz 
peek PKL SSI = il-oy e > T.) ui Ete UE or OL 
2) tk 
pORceae (eit ees) 


The parameter, zZ, is obtained from geometry considerations. 


For a cylindrical fiber shape, the expression for z is, 
2 
Z = 7d/2D ORT . 29> 


iemone—-hnalt factor in equation III.29 accounts for the un- 
certainty of the actual amount of exposed fiber surface area, 
which resulted from the initial combustion of epoxy. The 
non-dimensionalization of equation III.28 is presented in 
Appendix A. 
2. Internal Flow Heat Transfer Equation 

The energy balance on the internal flow is obtained 
by a procedure similar to the energy balance on graphite. 
Figure 5 presents a "smeared" differential volume of air, 
pdv, which shows the heat transfer mechanisms. The heat 


transfer mechanisms are: gq heat conduction through the 


cond’ 


air; q.¥4.: convection from the fibers to the air; mh, 
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energy transport due to the air flow, where h is the 
enthalpy of the air. 


Performing the energy balance (neglecting the higher 
order terms) yields, 


x ag 
z cond 
Ioond ah doonv Ioond ox dx 


dE 
aE pdxdydz 


Gite =s ) 
Cancelling terms and rearranging, the expression becomes, 
>Foond ° 3h 
-— ——— dx - m — dx 
ox a dx 


= og 
cea OE pdxdydz (oigake o3).) 
a. Thermal Conduction 


As before, Foo 


nq Was replaced by Fourier's law 
of heat conduction in the form, 


oT 
= = 8 
ee a Kk. xs dydz ieiele3 2 ) 
where Ke tomente  CcOMGUGCETVLtY OL air. 
b. Energy Transport by Flow 
Using the perfect gas assumption for air, the 
enthalpy h was expressed as, 
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where Cc. 1s the specific heat of air at constant pressure. 
The specific heat is treated as constant since it does not 
vary appreciably over the range of temperatures associated 
with the problem under investigation. It follows that, 


aT 
a 


= c, =3 (E38) 


Assuming the effects of the combustion by-products on air 


are negligible, continuity considerations give the mass 


flow as, 


m= Pa Up dyd2 (ticket, 34) 


where the density of air and the pore velocity are evaluated 


at ambient conditions. 
c. Convection Heat Transfer 


The convection heat transfer, — given by 


ealeatton I[11.23, with h; obtained from equation III.24, 


becomes 


oe eT) a (III.35) 


d. Change in Internal Energy 


For a perfect gas, the internal energy is, 
abe = poxezdT 
aa oa 
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where co is the specific heat of air at constant pressure. 


From this it follows that 


dE os 
== 2 dxdydz = PCA wee P dxdydz ELT. 36) 


Substituting equations III.32~-III.36 into equation 
III.31, the heat transfer equation for the internal flow 


becomes, 


5 oT. . oT. 
(Ke re pdxdydz = WS Ws. p dxdydz 
she 
- h; Nn = ie) = ec, 3E P dxdydz 


As in the energy equation for the fibers, the 
porosity, p, accounts for the fact that some of the volume of 
the porous medium is not air. Dividing the last equation 


through by pdxdydz and letting z denote dA/dxdydz, yields 


5 aT. ; oT. hz 
3x'*a Tx) 7 Mcp Ox * Tp Ve Be) 
oT. 
= eae 7t Clue 3 7 ) 


Non-dimensionalization of equation III.37 may be found in 


Appendix A. 
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S-enmyegen Transport Equation 


The final consideration in the formulation of the 
model is the mass transport of oxygen. The mass transport 
equation is obtained by a mass balance on a differential 
volume of porous medium. In accordance with convention, 


the mass balance is given by 


Oxygen into _ Oxygen out A Oxygen x Oxygen 
dv of dv Consumption Accumulation 


Figure 6 presents a differential volume, pdV, showing the 
relevant molecule-mass transport mechanisms. The mass 


transport mechanisms are: m, is the molecular diffusion; 


B 
mo is the convective transport due to mass flow; et (T 4) 
is the consumption of oxygen due to combustion; r is the 
reaction rate, and f is the stoichiometric ratio of the 
reaction. Representing the terms in Figure 6 by Taylor 


series expansions, the molecule-mass balance becomes (neg- 


lecting higher order terms), 


, om. om 
Mp + ma = Mp + 3x ax + mM . ae 7x dx 
i a 
i Fx (T,»¢) pdxdydz + =ppdxdydz Gini } 


Upon cancelling terms and rearranging, equation III.38 


becomes 
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3m 3m 


= se ed aos = po 
—— dx a bye Bee ney. C2 se PAxdydz Chie. 39) 
where » is the concentration of oxygen. The individual mass 


transfer terms in equation III.39 will now be developed. 
a. Molecular Diffusion 
Molecular diffusion m,, is given by Fick's law, 


mn = - B =o dy dz (TIT.40) 
where B is the diffusivity of oxygen into the porous medium. 

The diffusion of gases in porous media is limited 
by two factors. One is the collision of the gas molecules 
with other molecules, and the second is the collision of the 
gas molecules with the solid pore walls. This second phenom- 
enon is discussed by Bennett and Myers [17] and is referred 
to as Knudsen diffusion. To establish which behavior dom- 
inates, the mean free path of the gas (in this case, oxygen), 


must be calculated. Using the expression given by Treybal 


[18], 


w = (3.2u/P) (RT/2ng_M) */? 


where S. is the gravitational constant, the value obtained 
for mean free path, w, is smaller than 0.16. Thus, the 
collision of the oxygen molecules with those of air will 


restrict the diffusion process. On the contrary, if the pore 
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diameter is smaller than the mean free path, then the colli- 
Sion of the molecules with the solid walls will be the 
limiting factor. In the intermediate case, the mean free 
path and the pore diameter are of the same order of magni- 
tude. In this case, Bennett and Myers [17] proposed the 


following expression be used to obtain a value of diffusivity, 


I 
— = = + Cen. 41) 


= 
R g BY, 


R 
based on molecule-molecule collisions; and BS is the diffusivity 


where B_ is the resulting diffusivity; Bi is the diffusivity 


based on molecule-pore wall collisions. 

For the problem under investigation here, prelimi- 
nary calculations of the mean free path of the oxygen showed 
that diffusion is primarily of the inter-molecular collision 
type. With respect to this, and knowing that molecular 
diffusion is highly temperature dependent, the diffusivity, 


B* was obtained from an expression presented by Gilliland 





Peo"). 
pi/2 
a a ly73 1/3,72,,,-l aie 7 2 Zz 
B 435.7 > (Vo +Vo5 ) (M. + My5) cm/s 
(iny , 42) 


where P is the absolute pressure; = and Vo2 are the molecular 


volumes of air and oxygen respectively; M, and Mo are the 


molecular weights of air and oxygen, respectively, and B* 


n~ 


is the diffusivity. For this expression, T. is in degrees 
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Kelvin, and B* has the units of em“/s. As a result of 
tortuosity and porosity affecting the process, Denbigh and 
Turner [11] proposed that the effective molecular diffusivity 


for a porous medium is given by 


B =e OB ~/t < 
The porosity p accounts for the void ina differential cross- 
section. 
b. Convective Transport 
The convective term in equation III.39 is given 
by 
m = pu, odydz (III.43) 
The pore velocity u, is a temperature dependent variable, 
and is so treated. 
c. Consumption Rate 
The reaction rate term, EATS?) was discussed 
previously as equation III.26. The consumption rate of 
oxygen is determined by multiplying r by I1/f. 
Substituting equations III.40 and III.43 into 


equation III.39, the oxygen transport equation becomes, 


9 ,B*® 99 3 
Be 5) weeudydz - 3x (Up?) pdxdydz 


~ (i my 
fat a) Boxdyaz = 3, pdxdydz 
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Dividing both sides by pdxdydz, and letting B = B*/t, the 


expression becomes 


eee ee Se Se 
Sete) teeing) ~ (z)r(T,,9) = 5b. (III. 44) 


Expanding the second term of equation III.44, 


lS 
* 1-O- 
+ 

Q 
Q 
ho 
e- 


9 _ 
3x 4p?) = th 


and noting that ue depends on Dis from equation III.12, the 


expression yields 


5 34 du oT. 
ax(Upe) = YU, zy t Tr Ox a ee 


on ae is obtained from equation III.12. 
Substituting equation III.45 into equation III.44, 
the oxygen transport equation becomes, 


du. oT 


e 99 pa, _ (lh = 
a YD OX aT. aX ? ie Eite?) 7 


-) 
-) 


Q 
5S 
3S 


(B (III.46) 


“i 
@ 
ct 


In contrast with the heat transfer equations, mass balance 
does not depend on porosity. However, porosity does influence 
the oxygen concentration by appearing in the boundary condi- 
tions. Non-dimensionalization of equation III.46 is pre- 


sented in Appendix A. 
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IV. BOUNDARY CONDITIONS 


It was intended that boundary conditions be imposed 
which simulate the composite plate "burn" experiments of 
Fontenot [1]. However, for a one-dimensional model, the 
boundary conditions can only be approximated at best. The 
difficulty is associated with the existence of momentum and 
thermal boundary layers on the exterior surface of the por- 
ous plate. Complications also arise from the "blowing and 
suction" effects which are caused by the flow through the 
plate. 

The following boundary conditions provide a reasonable 


model for the one-dimensional problem, 


aul’ 


ao - _ nt m4 = 
(Ko +k) y hy(T,-T,) + oe(TZ- 7) at x= 0 (Iv. 
ot g 7 
Sey ae Er Fg Whe, = Bee = SE ke - re) ae x = & XIV. 
J8 eh at x = 0 (IV. 
a co 
aT oT 
———. eee at xe=L (IV. 
aX aX 
B a = u_(o - po_) at x = 0 (IV 
x p © : 
J = 
oe 0 elie ea EY 2 
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Boundary conditions IV.1l and IV.2 account for convection 

and radiation heat transfer from the porous solid. With 

an influx of air at x = 0, the heat transfer coefficient, 

hy, will depend on the magnitude of the flow entering the 
plate. Neglecting "blowing and suction" effects, several 
empirical expressions based on free convection were used to 
obtain a heat transfer coefficient at x = 0. These expressions 
yielded a range of values for hy from 1 tems Btu/hr.ft*-F. 

AS an approximation, hy was taken as the average value, 2.0 


2 


Btu/hr-ft":F. However, choosing h, as any value in the 


1 
above range did not affect the results of the analysis. 
The magnitude of the forced convection heat transfer 


coefficient at x = L, h varies in a direction parallel to 


Ps 


that of the external flow, U_. In addition, h, depends on 


2 
the efflux of the gas at the surface. To simplify the analy- 
Sis, h. waS approximated by the relation for a smooth flat 


plate given by Holman [20] as 
(avi. 7) 
for laminar flow, where Pr is the Prandtl number and L* is 


a reference length, arbitrarily taken as unity. For turbulent 


flow, 


nw 


8 


Pee (oe 7 Res — 850) . (IV.8) 


Nd 
Es 
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Kays [21] provides an alternate scheme for treating boundary 
layers and "blowing and suction" effects. 

Boundary condition IV.3 accounts for the air entering 
the porous plate at ambient temperature. The Dankwerts 
boundary condition proposed by Riaz [6] which accounts for 
conduction-convection interaction, is a more reasonable 
simulation. 

The boundary condition at x = L for the air heat trans- 
fer equation appears to be a reasonable one. It follows the 
behavior observed in the interior of the plate and does not 
fix the temperature of the air at the wall. The boundary 
eendition, aT / 9x = 0, used initially provided results simi- 
lar to those obtained using equation IV.4. 

Boundary conditions IV.5 and IV.6 are Dankwert conditions 
[22] for flow through porous media. A convincing discussion 
for the Dankwerts conditions is given by Bischoff [23]. A 


brief summary of the discussion is presented in Appendix D. 
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V. FINITE ELEMENT METHOD 


A. GALERKIN FORMULATION 


A Galerkin formulation of the Finite Element Method was 


used to obtain solutions of the field equations. A conven- 
ient form of the field equations was used in the formulation. 
These equations are presented in Appendix A. 
@eephite energy equation 
00 
3 30 7 = —  (V.1) 
any rw, Ao (0 Tr) + api seat) = 4 Pe) 
air energy equation 
3 ar ar 7 ar 
Fal on ae an + v3 (0 rT) oS | at (We? } 
oxygen diffusion equation 
3 3 > 30 ait " 30 
wer an oT es on wi ~ 5 dt a) 
where the coefficients .are defined in Appendix A. 
Field equations III.28, III.37 and III.46 subject to 
boundary conditions, IV.1 - IV.6, and initial conditions 


define the problem. 


The closed domain (0,L) was partitioned 


into (n-1l) contiguous elements of variable length Ray ee 1), 


--.-,(n-l). This defines an n nodal point model. 


SL 


The three 





non-dimensional field variables 0, [, % were approximated 


by 
O(n,t) = ¥,(n,t) = 2G, (n)®, (t) (V.4) 
Pim,t) = ¥j(n,t) = 2G, (n)T, (t) (V.5) 
o(n,t) = ¥,(n,t) = JG, (n)¢, (t) (V.6) 
where G., for i=l, ..., n, is a set of specified basis 
functions with local support, and the sets {O,,T,,%7i=l, eke; 


n} are the solution coefficients to be determined. The G. 

were selected to satisfy the condition Gj (ns) = AGG where 

the kronecker delta 845 is defined by ane = 1fori=j, 

and re = 0 for i#j. As a result, 0, Tf, and 6 are the 

values Yar Yor ¥ 3 at the nodal points (l.e., 0; (t) = ¥, (n3,t))- 
Linear interpolation functions (Figure 7) were used 

as the basis functions. These are the lowest polynomial 

functions which provide the necessary function continuity. 


As a measure of error, a residual function, rs, is 


defined for each field equation by 
r;(n,t) = A, (¥) a 2 = M2 3 iy ..7) 


where A. denotes the spatial operator of the ith equation. 


For convenience, the following convention for differentiation 


a2 





is adopted, 








For field equations V.1, V.2, and V.3 the residuals are: 


n n 
1=1 i=l 
nr e 
+ Lar (00) aie ‘I Go; (v.8) 
i=l 
n n n 
= t ! rae t = 
1=L1 1=1 1=1 
n : 
- v4 (tT) oD GT; (V.9) 
1=1 
n n n 
Ys [w, (1) ) Gio,]' - wolf) 2 Gio, —w3(T) 2 G.T.G.o5 
1=L1 i=l 1=1 
n r 
- w,r(0,%) - we wey Go, (V.10) 
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where the coefficients multiplying the response variables 

are themselves functions of the response variables, and thus 
the equations are nonlinear. In accordance with the Galerkin 
method, the final system of ordinary differential equations 
was obtained by setting each residual, Bac orthogonal to 


each basis function, Gs, that is 


1 
jeeceercne — 90) 1 = 1,2,...,n; (7nd) 

Ve? 
7 le a 


The 3n ordinary differential equations given by equations 

V.1ll retain the character of the original set of partial 
differential equations. Thus linear field equations trans- 
form to matrix operators and nonlinear, coupled field operators 
become nonlinear, coupled algebraic operators. Incorporation 
of the boundary conditions resulted in 3n nonlinear, coupled 


ordinary differential equations 
A(t) ¥(0,T,¢) + F(t) = B(t) ¥ cee) 


subject to initial conditions, where B is a 3nx 3n matrix, 
A is the operator associated with the field operator A; in 


expression V.7, and F(t) is an excitation vector. 


Adopting the convention, 
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and applying the operation of expression V.11 with an inte- 


gration by parts on the second order derivatives gives, 


il 1 
{G,},<G, >"{o} - ry geo! 


il il 
- pe apie’ | um ie ey (V.13) 
q a _ il 
az 3 ‘ {G; }dnr (0,9) ee hy f NG BSG ites 
0 0 
as 1 
(Gy}vy<Gartir}] - Vy f {G,}'<G,>'dnir} 


0 0 


1 1 
zsh Vo ooh Gade t + V3 gma sCarcn ie} (V.14) 
- - ue 1 ; 
Oe * ae ee SU et) 
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’ . 
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eG uw. <G.>'{ >} - w 
ah aoe Seats 


ie 
1 f {Gy} '<G5> tanto 
0 


1 il 
- W, fei (9) - Wa Mee TG iret) 
2 il _ 1 
= Poe tae aiee = Us pei e yr enio) (v.15) 


The first term in each of the above expressions is a boundary 
term which permits the incorporation of natural boundary 
conditions which will be shown in Section V.B. Coefficients 
X, V, W™ are variable-dependent properties, and were taken 
as the average value of the properties over an element. In 
the limit, as the elements get smaller (i.e., n+ ~), the 
average values converge to the exact values. 

Inspection of expression V.13, V.14, V.15 shows the five 


operators, 


It 


per <8? *én (V.16) 
1 
a ie (v.17) 
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1 
OP aides (V.18) 


1 
Se 2G, > dr (V.19) 
1 

Ps (V.20) 


To formulate these operators, the global shape function, G., 


was defined on the local level by 


= ic 1) (1) 
Gs a =)5) ® 39 vere 1) 


where Sy and G5 were defined by 


(1 - =) for € in element (e) 
e 
(2) 


0 for — not in element (e) 


for — in element (e) 


o | 


a. 


0 for &€ not in element (e) 


oe 





and de is the length of the eth element. The ©@ notation 


in expression V.2l1 means that G, Eomee UnLon | OF oe 


The local (element) shape functions have the 


and jae 


following properties, 


My ; 
GQ) f OP He = 0) if Gem 
0 
u aL Sr l= jj 
o (e) = - 
(il) g, (3) = jne = 


Having defined the local shape functions, the elemental 
matrix operators contributing to the global matrix operators 


V.16 through V.20 are, 


il -1 
1 
a 
ji ie2 <G5>'dr » i (V.16') 
-1 1 
-1 1 
1 
i] ie v 
J {G3 }<G,> dn B® 5 Cee) 
| -1 1 
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1 be 
f iG, }<Gj>dn > = cele) 
1 2 
ae =, 
ee ee 
1 Oe ak. 2 
-1 ' 
e226, > }<G5>a0 >» = id (V.19") 
i-l 1 
i. 
1 Le 
J (G tan > = (V.20') 
1 


The derivations of these operators are shown in Appendix E. 

Since there is extensive coupling in the field equations 
and there are three degrees of freedom at each nodal point, 
the numbering scheme for the system matrices was important. 
To minimize the bandwidth of the matrices, the numbering 
scheme represented in Figure 8 was used. 

Figure 9 represents the matrix A(t) of expression V.12 
for any three successive nodal points. The distribution of 
the elemental matrices for the FEM operators is shown, as 
well as the possible locations for the coupling of the field 
equations over an element. In addition, Figure 9 shows the 
bandwidth that would be observed for any n-1l element solu- 
tion. For this scheme, the bandwidth is nine. If coupling 
is not present, the bandwidth is seven. Figure 9 reflects 


the extensive coupling that is present in the model. All the 
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matrix elemental locations are filled except those which 


are shaded in. 


B. IMPLEMENTATION OF BOUNDARY CONDITIONS 
Having formulated the system matrices for the field 
equations, treatment of the boundary conditions will now be 
discussed. Each field equation is considered individually. 
1. Fiber Heat Transfer Equation 
The fiber heat transfer boundary conditions in non- 


dimensional form from Appendix A are 








5 oa - oeL 4 74 
n = 0 hy ae h,L(0 1) + F(T, ce) (Vie2 2) 
7 O— = _ oeL,-4_ 74 
nel hy in hoL (0 1) ep (Ty ae) (V.23) 
Since the first term in expression V.13 is 
iL 
t 
Mee Ws Sere {o} : 
or in analogous form 
iI 
EAC, 
dy ae ; (Va ) 





natural boundary conditions V.22 and V.23 may be directly 


substituted in equation V.24. 
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1” pet Ty! 


changing with time, are evaluated at the previous time 


The response dependent parameters, h h 


step. Thus, the boundary conditions are incorporated in the 


system matrices as follows, 


(ls) “hi: added to the stiffness matrix A(t) 


at location Ari 


(2) hyL - Sarg - T°): added to the excitation 


vector Gee at location Fi 





(3) -hjL: addea to the stiffness matrix A(t) at 
LoecaElon A3n-2,3n-2 
(4) hot - SS2(7t - 72): added to the excitation 


tT OCO@ 


vector F(t) at location F,__, 


2. Internal Flow Heat Transfer Equation 
The non-dimensional boundary conditions presented in 


Appendix A for the air energy equation are 


n = 0 ros 6] (Viz ) 
: eT 

n = | 7 - (V.26) 
The essential boundary condition at n = 0, r(n = 0) = 1 is 
imposed in the Galerkin equation as follows. The A, , row 

a 
of the A(t) matrix, the Bo ,; row and the B, , column of the 
~ f f 


B matrix, and the EF. location of the excitation vector, F(t), 
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are all set equal to zero. The Bo 2 location of the B 

Matrix is then set equal to one. The natural boundary con- 
dition at n = 1 for the air heat transfer equation was treated 
in the same manner as that of the fiber heat transfer equation. 
Foregoing the individual steps the boundary condition is 


implemented by 


hobv, 
(1) - ——— added femthe stittness matrixyA(t) 
1 ~ 
cee tocation een? 
ies a) 2 .: 
2) 21 _ i (p4 — T*): added to the excitation 
ST eT g <9 


vector F(t) at location F.- 


Si. Oxygen Transport Equation 


For the oxygen diffusion equation, the boundary 


conditions in non-dimensional form from Appendix A are: 

a . 

an Wo (¢ p) vee? ) 
= ite 

n aE ee 0 (V.28) 


Since these are both natural boundary conditions, they were 
substituted for the first term in expression V.15 at n = 0 
and n = l, respectively. Coefficients Wy and WwW, are evalu- 


ated at the previous time step. The boundary conditions are 


implemented by adding 
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yee towne SCIiTINess matrix A(t) at location 


~ 


ae 
and 
(P2) WoP: to the excitation vector F(t) at location 


— 


This concludes the discussion for the implementation 
of the boundary conditions. A word of caution is in order. 
After each time step integration, the time-dependent coeffi- 
cients of the boundary conditions (i.e., hy. ho, Kop. feats) 
are reevaluated. Before incorporating the updated coefficients 
iieomene Stiffness matrix A(t), the previous values must be 
subtracted out. The results of not taking this into account 
will be obvious. 

The implicit system of ordinary differential equa- 
tions was integrated numerically by a modified implicit-Gear 
method developed by Franke [2]. The time-dependent coeffi- 
cients of the FEM operators in expression V.13, V.14, and 
V.15 were updated at the previous time. The reaction rate 
term (III.26) was also evaluated at the previous time, and 
M@erears 2m the excitation vector F(t). In this way, the 


final system of ordinary differential equations was 
A(t) ¥ + F(t) = B i 


where A and B are (3nx 3n) matrices of temperature dependent 


~ 
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coefficients. F is a (3nx1) vector arising partly from the 
reaction rate terms and partly from the boundary conditions. 
Bemienea, the elements of f are dependent on both temperature 


and oxygen concentration. 
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VI. RESULTS AND CONCLUSION 


A. NUMERICAL CONSIDERATIONS 

The input parameters which were varied in the computer 
analysis include (1) plate thickness, L, (2) fiber diameter, 
d, and (3) air flow rate, U_. Other parameters which could 
be varied are the ply thickness, D, ambient temperature, T , 
and fiber emissivity, «. One set of initial conditions was 
used for the analysis. These are shown in Figures 10 and ll. 
For this initial effort, the actual initial conditions were 
not of prime consideration. The selection of initial condi- 
tions of Figures 10 and ll will be discussed in Section VI.C. 
The finite element program calculates the remaining system 
parameters such as permeability, porosity, pressure differ- 
ential, and the response variables as functions of time and 
position. Parameters which are functions of temperature, 
such as Kos Kg, Kos h., Os: u and yp are continuously updated 
during the transient analysis. 

The preliminary solution effort showed the system of 
equations to be very stiff (refer to Shampine/Gordon [24] 
and Gear [25] for a discussion of "sStiffness"). Changing 
the integration algorithm from a sixth-order Runge-Kutta 
(IMSL subroutine DVERK) to a modified implicit-Gear method 
for stiff systems, developed by Franke [2], resulted in a 


Significant reduction in CPU time. The computational effort 
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was carried out on an IBM 360/70. Typical runs required 20- 
30 minutes CPU time for problems for which ignition occurred, 
and 6-25 minutes CPU time for extinction problems. Ignition 
problems were terminated when the temperature at any nodal 
point exceeded 2200 degrees Fahrenheit, or when the graphite 
at a nodal point was totally consumed. Extinction problems 

. were carried out to steady state. 

A twenty five nodal point model (75 o. d. e.) provided 
results which differed less than five percent from the re- 
sults of a 32 nodal point model (96 o. d. e.) and was adopted 
for all computer runs. For the twenty five nodal point 
model, approximately 275K bytes of core was required. The 
computer program which includes the FEM formulation and the 
integration routine as well as a sample input file is pre- 
sented at the end of the appendices. 

The numerical solution produced satisfactory results with 
one exception which will now be discussed. For a particular 
range of temperature and oxygen concentration (approximately 
1500 degrees Fahrenheit and .001 lbm/ft”, respectively), 
there is a significant increase in reaction rate. This 
accelerated reaction rate produced a negative/positive 
oscillation for the nodal values of oxygen concentration. 

The oscillation occurred in a region of the plate where the 
oxygen appeared to be totally consumed. As discussed by 
Frank-Kamenetskii [3], nth order reactions may yield zero 


values for the oxygen concentration and gradient within the 
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medium. The cause of the oscillation was as follows. 
Although the reaction rate is a function of temperature and 
oxygen concentration, numerically, it is treated as constant 
during an integration interval. As a result, the concentra- 
tion in the region of high reaction rate becomes negative. 
The reaction rate was updated for the next time interval by 
setting negative concentrations to zero. During this time 
interval, for the region of zero concentration, there is no 
reaction. Without consumption, the oxygen concentration will 
increase and the cycle repeats itself. The oscillation was 
aggravated by large plate thickness and low permeability. 

The instability occurring for the high reaction rate may be 
corrected by (1) numerically integrating the rate term in 

the time domain, (2) iterating until convergence is obtained 
between consecutive values of concentration, or (3) decreasing 
the time step and the length of the elements. As a matter 

of expediency, measure (3) was used to minimize the insta- 
bility. The fiber and air temperatures were essentially 


unaffected by the oscillations in the oxygen concentration. 


B. RESULTS AND OBSERVATIONS 

Table 1 presents the results of fourteen problems. In 
Table 1, the parameters which depend upon temperature a? 
Re,, and h;) are evaluated at ambient conditions for com- 
parison purposes only. These parameters, in fact, will vary 


during the course of the transient analysis. In all cases, 
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the initial condition on graphite temperature was 1050 
degrees Fahrenheit. 

Figures 10, ll, and 12 show the transient behavior of 
Case l. For the initial conditions shown, ignition occurred 
in approximately eight seconds. The reaction rate and oxygen 
concentration responses for Case l are typical of the problems 
for which ignition resulted. Figure 13 shows that the air 
and fiber temperatures for Case 1 do not differ significantly. 
This behavior is typical for most problems. However, for 
plates with high porosity (greater than .8) and subjected to 
relatively high flow rates, significant differences in the 
air and fiber temperatures do occur (see Figure 14 for Case 
9). In a previous effort by Vatikiotis and Salinas [26], 
linearly varying initial conditions were investigated. Although 
the reaction rate used in that investigation differed from 
the one used in the present investigation, the overall re-~ 
sults remain valid. The main observation was that ignition 
is less likely to occur for thin plates with the ambient air 
entering the hotter plate surface. Since temperature con- 
trols the reaction in the kinetic regime, cooler air entering 
the hotter plate surface enhances the heat loss. This reduces 
the fiber temperature, thereby decreasing the reaction rate. 
In effect, the cooler air enters the plate where it is most 
needed. 

In the present analysis, observations were made by vary- 


iiesonesinput spearameter (1.e., Lor d or U_) while keeping 


68 





all others fixed. The effects of each parameter on the 
behavior will be discussed individually. 
1. Effects of Exterior Velocity 

As shown in Figures 15, 16 and 17 (cases 4, 1 and 
2, respectively) for an initial condition of 1050 degrees 
Fahrenheit, three distinct regions of combustion were ob- 
served (extinction-ignition-extinction). These became appar- 
ent as U_ was increased from low velocities (10 knots) to 
higher velocities (120 knots). Vulis [4] discusses an experi- 
ment of a heated carbon rod with an air jet impinging upon 
it. For a certain range of flow velocities, ignition was 
observed. However, extinction did occur for velocities less 
than and greater than the velocity for which ignition occurred. 
The behavior of the carbon rod and that of the graphite mat 
is similar. This behavior may be explained by the Semenov 
Hegel of Figure 1. For high U,, the internal heat transfer 
coefficient, hi, will be large as a result of an increase 
in the pore velocity. Since the slopes of the heat loss 
lines, Gas increase with increasing h; for a constant To 
the graphite temperature decreases from the critical ignition 
point I. This causes extinction for the higher exterior flow 
velocities, U,. In contrast, for low U_, the slope of or will 
be small and ignition is more likely to occur. However, the 


heat generation curve, om changes position due to a decrease 


in oxygen concentration. This reduction in oxygen concentra- 


tion is caused by the lower induced flow rate through the 
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plate. The overall shape of dg is flattened such that the 
qo line crosses Ss behind the new critical point I. Thus 
extinction is observed for this region of low U_. 
2. Effects of Fiber Diameter 

Varying the fiber diameter with the remaining input 
parameters fixed, also yields three regions of combustion 
(extinction-ignition-extinction) for a fixed initial tempera- 
ture. The results are shown in Figures 16, 18 and 19 (Cases 
1, 9 and 10, respectively). The fiber diameter and ply 
thickness determine the permeability of the plate. Similar 
to U_, permeability affects the pore velocity (i.e., low 
permeability + low velocity). In turn, this affects the 
convective heat transfer and the amount of oxygen entering 
the plate. These three regions of combustion are explained 
by the Semenov model for the same reasons as those of varying 
Oe 

Porosity 1S a measure of void space and denotes the 
Space available for oxygen. Porosity 1S associated with the 
internal geometry of the medium, as iS permeability. However, 
they are basically different parameters. Permeability 
(hydraulic conductivity) is associated with convective air 
flow through the medium. Thus, an n-fold increase in the 
Fiber diameter and ply thickness will not affect a change 
in porosity, but will affect an n? order change in the permea- 
bility. For ann greater than one, the result is a decrease 


in pore velocity without a change in porosity. Porosity 
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determines the maximum oxygen concentration per unit volume 
of fibrous mat. This indirectly affects the reaction rate. 
Thus, for a given u,s the lower the porosity, the less oxygen 
there is for combustion. It appears that high porosity 
would enhance combustion by the presence of more oxygen. 
However, higher porosities also provide more fluid per unit 
volume resulting in greater heat transfer. This is observed 
in Case 9, Figure 18. 

3. Effects of Plate Thickness, L 

For the given initial temperature of 1050 degrees 

Fahrenheit, varying the plate thickness yielded two regions 
of combustion (extinction-ignition). Extinction was observed 
for thin plates (Case 2, Figure 17). As shown by Figures 20 
and 21 (Cases 7 and ll, respectively), ignition was observed 
when plate thickness was increased. Plate thickness also 
affects the pore velocity. Decreasing L, increases the 
pressure gradient across the plate, thus increasing the pore 
velocity. The result is that extinction is more likely for 
thin plates because of the enhanced convection heat trans- 
fer resulting from the increase in pore velocity. Whereas 
pore velocity is inversely proportional to L, it is an exponen- 
tial function of U_. For the range of parameters in this 
Mivecetgation (1.e., .1” < L < 3.", and 10 knots < U_ < 100 
knots), an n-fold change in U_ will produce a greater change 
in ee then will an n-fold change in L. Compare Cases 1 and 


2 to Cases 7 and 8. 
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PeGcbearmacteristic feature of ignition for thick 
plates was observed. As the thickness increased, the igni- 
tion region (the spatial location) moved closer to the 
entrance surface of the ambient air. Compare Figures 20 


and 21. 


Gee D> LSoCUSSION 

At this point, discussion is made of ignition tempera- 
tures for porous graphite plates. The ignition tempera- 
ture for Case 1 was 1050 degrees Fahreheit. This tempera- 
ture was obtained by varying the initial temperature until 
ignition occurred. In all subsequent cases, 1050 degrees 
Fahrenheit was adopted as the initial condition, and we 
observed whether extinction or ignition occurred. The igni- 
tion temperatures given here are the result of taking n = 5 
in the nth order reaction rate. Significantly different 
ignition temperatures will result for other n. For example, 
for n= 1 (lst order reaction), the ignition temperature for 
Case l was approximately 1600 degrees Fahrenheit. As Frank- 
Kamenetskii [3] has suggested values of n between 1/3 and 
2/3, the average value of 1/2 was adopted for this analysis. 
The ignition temperatures obtained using n = > agree favorably 
with the experiments of Fontenot [1]. For fiber combustion, 
ignition temperatures for porous graphite plates are on the 
order of 1100 degrees Fahrenheit. 


Consider the results of varying U_. Extinction occurred 


for Case 2 and Case 4; thus, the ignition temperature 
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for each case is greater than 1050 degrees Fahrenheit. 

Since ignition does occur at 1050 degrees Fahrenheit for 
mieermediate values of U_ (Case 1, Case 7), ignition tem- 
Pemarure is not a monotonic function of U_. The precise 
determination of ignition temperature for a specific value 
of air velocity can be obtained by a large number of com- 
puter runs. A similar discussion can be made for the igni- 
tion temperatures dependence on fiber diameter. The general 
shape of these functions is represented by the shaded region 
on Figure 22. These suppositions may only be valid in the 
limited range of parameters associated with the present 
investigation. 

It is interesting to observe that changing the initial 
condition of the oxygen concentration only affects the early 
part of the problem. The response of the concentration is 
fast compared to that of temperature. After approximately 
one second, the transient behavior is the same regardless 
of the initial condition on the oxygen concentration. 

A brief description of the Semenov model was given in 
Section II.A. That discussion was for a lumped parameter 
model. The actual combustion process is more complex since 
there are spatial variations in temperature and in oxygen 
concentration. The results of this analysis show that both 
kinetic and diffusion regimes may exist simultaneously in a 


porous medium. Hence, it is not reasonable to restrict an 
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analysis of combustion in a porous medium to the kinetic 


regime. 


D. CONCLUDING REMARKS 

The model has provided considerable insight as to the 
behavior of fibrous composites subject to combustion. 
Improvements may be realized by extending the model to 
include: 

meee tne: Combustion of the epoxy Matrix (1i.e., start 

problem with first stage of combustion process). 

(2) The change in geometry of the porous medium due 

to combustion (affects p, m and 2). 

(3) A more complex reaction (i.e., secondary combustion 

of gaseous by-products, such as carbon monoxide). 

(4) The effect of gaseous by-products on fluid properties. 

(5) Generalization to two and three dimensional models 

(results in anisotropic properties). 

The results show that certain properties, such as plate 
thickness, permeability and porosity, have significant effects 
on combustion. Therefore, one could select these design 
parameters to improve the flame resistance of composite 
Materials. This would be beneficial to the survivability 
of a composite structure. A comprehensive set of analyses 
could be undertaken to couple the combustion problem to the 
strength problem. In this way, one could achieve a design 


which maximizes strength and minimizes combustion. 
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Another use of the present model would be in the study 
of energy systems utilizing particulate fuels (i.e., coal, 
biomass). In this case, the design parameters could be 


selected to maximize the performance of the system. 
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APPENDIX A 


NONDIMENSIONALIZATION OF FIELD EQUATIONS 
The field equations and boundary conditions are: 


aT Rez oT 
ca | Sy . 2 Z BA _ J 
wal! (kotk,) oe (os oe T,) “+ Aur (T 7%) oe aE (Lene 28) 


5 oT. ; oT hz oan 
me aex. “Saox pg ‘a’ = “aca ve (ITT. 37) 
Ju. oT 
9 net) ~ y 9% ~ _Py__8), - hb £23o 
ax (BSy) 7 UDsx aT. ‘Gx’? (FETS 9) = FE (IIT. 46) 
at x = Q 
of ~4 74 
cise oe = ee + SES ee) Civ 1) 
= T. is) 
B ey = u_(¢- ) (iS) 
9x = Dp ) po. ° 
at x = L 
oe ~4 24 
(Kot) ae Sieh NE ee - em (IV.2) 
aT. aT 
ox =~ (IV. 4) 
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The 


Q7} @ 
 |-o 
! 
© 


LIV. 6) 


| 


The nondimensional variables are defined as: 


6 = Ba Be nondimensional fiber temperature 

rT = T/T, nondimensional air temperature 

d = $/¢, nondimensional oxygen concentration 
n = xX/L nondimensional distance 


time variable, t, will not be nondimensionalized. 


Using the temperature of the fiber, Tye to demonstrate 


the technique of transformation, we have, 


Al = T 0 
Gg co 
It follows that, 
aX o OX 
3 oT 2 
sane Se = hy 3 90 
2 oo 
ox OX 


Using the chain rule to transform x to n, gives 


ong = op 20. an 
ax o dn @X 
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Since n = x/L, the partial derivatives become 


org = #20 
aX L n 
and 
3° Tr T .2 
Moe 8 on 
L 20x 
aX an 
Therefore, 
3° Tr T 2 
ee 8 
2 ly) 2 
aX L an 


Similar transformations for T. and ¢ yield, 





“Se ee ee ag _ fa ae 
ox 1 2 elag ax L an 
sn T 2 Ae. 
2h = pees a eo *. 
3x? pe on 3x7 L* an 


Substituting these relations into the field equations and 


removing the constants from the differential operators gives, 





0 3 30 i 0 
ae oe ae _ = = a 
ae ant {Kg LS ine. Tiss” ee 2 eat Pg gat 
Mee eet, i.) - , co ot 
r2 an aan a on Dp aa vat 
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d eee ¢ T du 
Memory PS Oe | ee, py er, _ (bh = ae 
L an Pan? L an L ST) an° ar) ‘ee dt 
Upon rearrangement, the equations become, 
lel 2L° 2 
fe) 00 1 AHL 2900 
ain | eee on. laepy Tr) + pe (Tyo) oar yt 
2 
glen, ral? 
: my (Ow i _ 2or 
mo aen) Catan * (ONS 2 oS ig 
du 2 
ee 8 ea (Byte - (hye - z2 22 
ieee pean foam an? ~ SFG 7 "Tg? = at 
Letting, 
2 
bine 
a Se: ets 2 
ee Sot Xy) oP ny @ i -) a a Tice 
h. zu 
v = k v, = me_L TS v ome ine 
Ji a 2 a 3 Pp 4 aa 
du 2 
= z = __P iL’ 
Wy B Wo Se W 9 TV ST. Wy ra é_ 
W = 12 
the field equations may be written as, 
9 96 - 00 
aT (Ay aa A, (0 T) + Agr (T +o) Sere: re) 
3 ar ar _ Oe 
ame ey, Vo am ate v3 (0-1) = A AT (A.2) 


io 





a: x: ar | 2 ad 
Se ee oer ston, aT Tg?) = 85 oe TAS 
Noting that the coefficients \’, v, w are response dependent 
variables, equations A.l, A.2, and A.3 are considered in 
their final form in Section III. 

Substituting the nondimensional variables into the 


boundary conditions gives, 








ea ” oeL 4 74 = 
Ay Pini = Ah Lo 1) + T a, tT.) x = 0 
. = aL ee 6 

ad i 7 
Wy er = wo (e p) x = Q 

and 

a0 oeL 4 “4 7 
hy a = hoL(o-1) - T_ ie = i) x = L 
a. 30 
an on “E 
0d _ = 
oh = Q x = L 


The dimensional form of the absolute temperature will be 


retained for convenience. 
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APPENDIX B 


POLYNOMIAL APPROXIMATIONS OF THERMAL PROPERTIES 


Relations giving the viscosity and thermal conductivities 
of air for varying temperature were required. A simple way 
to obtain values for these properties is to fit empirical 
data with 2nd order Lagrange polynomials. 

The general form of the 2nd order Lagrange polynomials 


for thermal conductivity and viscosity are 


(T7Th2) (T-T, 3) ee ae ela! 


k = k ———ooSoO k 
(7) th ae) a = 
- ToT Ta2) Tar 733) al (Too T43) ‘Ta Tal? a2 


ee ee: 


ili 

te i = Lew oh a k (B.1) 

73 ta? ‘743 T 32? a3 
—_ (T.-T,2) (T,-T,3) ae (T.-T,3) (T6-Ta1? ' 
z (T..=T .)(T .-T .) Vitae.) nT =|) 
(To Ta2) Tay tae 1 (T2 T43) ‘Ta Ta? 2 

Tata?) Tan TMa2? 

* [€_,-T_,)(T_.,-T_,) 43 ooo) 

aor al a3 a2 


Choosing three points from a range of temperatures that 
would be representative of those observed during the analysis, 


the corresponding values of the properties are: 
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(subscript) 


Temp k u 


a 
Btu lbm 
(F) “Tec earee Esa 
0 en 0394 
700 0284 0765 
1500 0423 107 


Substituting these values into expressions B.1 and B.2 


gives, 


a 


u 


eavuseSisnnnes °°) © + @is00=0)11500-700) 


(T-700) (T,-1500) CL ~1500) (T_-0) 


(0-700) (0-1500) ee) * FOG UOC) (HOO) 


(To -0) (TL = O10!) 


+ TreqocoyTIs00=700) (0423) 


(T_-700) (T_-1500) (T -1500) (T_-0) 
a a e = (.0765) 


Cr 4 0) (TL 7 O10) 


+ reqs) is0d=7007 (107) 


These expressions reduce to, 


k = 
a 


u — 


SO) = 5 Btu 
-1.041x1072r? + 6.029x 10-7 + .0394 LOM 
a a ft-hr 
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Both expressions give property values which are within 
five percent of the data for temperatures to 2000 degrees 


= 


Fahrenheit. 
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APPENDIX C 


TRANSFORMATION OF THE REACTION RATE TERM 


The expression for the oxidation rate of graphite fibers 


taken from Parker and Hottel [15] is 


6 
— = Pas exp (744000) = (C.1) 
J T R T cm ~sec 
gk gk 


where R is the universal gas constant (1.958 cal/gmole-k) 
and Poo is the partial pressure of oxygen (atm.). The 
partial pressure may be represented in terms of the tem- 


perature, T, and oxygen concentration, ¢, by the Ideal Gas 


law in the form given by Kanury [16] 


M x 


a 
P = (——) R, T 9 (CE) 
02 Moo A 
where Mi/Mo> = .9 is the ratio of the molecular weights of 


air and oxygen; R, (= 53.34 ft.-lbf/lbm R) is the gas 


A 
constant for air. Substituting Mi/My > into equation C.2 


and converting P ) to atmospheric units the expression 


0 


becomes 


p A See Ree? G atm. 
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Converting es to lbm/f£t*-hr from gm/cm“-sec, equation C.l 


j 


becomes 
LO 
E _ 7.04 x 10 p Soe) lbm 


euestituting in for P and converting the absolute tem- 


02 


peratures from Kelvin to Rankine, the rate expression 


becomes, 
R, T $ 
r= 4.01410 5 exp (—S***) = 
g T Je ie ee ga 


Z| SI 


The absolute temperature at the fiber surfaces is Bee Thus, 
T in the numerator becomes Ty! and the expression becomes 
7 2 


ig = 4.014x10 R, T 
S) A g 


pe) cree. 


st 
g 


d exp ( ) 
To transform the reaction rate to Tyagi y, the ratio 
of a fibers surface area to its volume is placed in the 


numerator of Ig: The expression for y is, 


= 
Qu 


(ex 3) 





2 
d 


KG 
it 
NS) fee 
>! = 
Ay 
i] 


where dis the fiber diameter, and the = factor accounts for 


Nye 


the uncertainty of the actual amount of fiber surface after 


the epoxy has burned away. The reaction rate of the graphite 
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fibers per unit volume is expressed as 


r= 4.014x10’ y R, Te’? 6 exp (284) Lb 
2 Ty ‘a laa 8 


To obtain the heat liberated per lbm of fiber, r is 
multiplied by AH, the enthalpy of formation of carbon dioxide 
per unit mass of graphite consumed. This value is approxi- 
mately 14,085 BTU per pound-mass of fiber and was obtained 


from Kanury [16] 
Heat generation rate = AH E(T +) 


To obtain the oxygen consumption rate, r is multiplied 


by the inverse of the fuel-air stoichiometric ratio, 


Oxygen consumption rate = (=) BA go) 
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APPENDIX D 


JUSTIFICATION OF THE DANKWERTS' BOUNDARY CONDITIONS 


The Dankwerts' boundary conditions for the oxygen diffu- 


sion equation are: 


qx =) UB} 7 PO.) 2S SU 
do _ 2. 
as 0 x = L 


Bischoff [23] presents a discussion of these equations. A 
brief summary is provided for completeness. 
Figure 10 shows the one-dimensional oxygen transport 


region considered in this analysis. 


x = 0 x= L 
Entrance Porous Exit 
Section Plate Section 
Flow 
ane? ae a te a Wen 7G 


The concentration of oxygen, ¢, for each region is distin- 
guished by subscripts; B and u are the diffusivity and the 
velocity, respectively, for each region. The oxygen diffu- 


Sion equations for each section are as follows, 


87 











B, d 

Pella « 4 x < 0 

uy dx? dx _ 
Bas _ do _ in ) = 0 Ca Se eal, 
u 2 dx ro pe a 

by dx 
B 2 d 
~2 dg _ “2 = 0 Dee aE, 
Us gy dx — 


with the general boundary conditions, 


$4 (-@) = >. 
pd, (0°) = 9(07) 
B, do, (0 ) S 
eee mee 7 i aeeeechon(O. 
Po, (0 ) Po ax SS u ax 
1 p 
po(L) = $5(L") 
Been i(t)” ) 
~p, B de(L) _ 0 Ce 
oo) P u,, ax - $5 (1 ) Us aa. 
do (+2) = vealialee 
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aye 


CDs 


UDG 


ot) 


2) 


3) 


Ay) 


72) 


.6) 


7) 


8) 


se) 





For convenience, the parameters U;, B; will be treated as 
constant. In addition, the porosity, p, has been included 
to account for the discontinuity that exist for the concen- 
tration. The discontinuity is caused by the oxygen from a 
totally gaseous environment to one that is partly occupied 
by solid. 
An analytical closed-form solution for this set of 

equations is not possible because of the nonlinearity of 
equation D.2. However, the solutions of equations D.1 and 


D.3 are, 


4 = A. +A exp (u,x/B,) Ce 0)) 


exp (u,x/B,) (De 11) 


5 


Applying boundary conditions D.4 and D.9, the following re- 


sults are obtained, 


The solutions D.10 and D.11 become 


= 6, + A, exp(u,x/B,) UB) Le) 


7 2 
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—_— (D.13) 


It would be necessary to have the solution for the non- 
linear equation D.2 to solve for A, and A.- However, the 
constants need not be known to continue with the analysis. 


From equation D.12, 


,(0 ) = $, + Ay 
and, 
do, (0 ) 7 a 
dx By 2 


Be ul + 
» lige: 2 + B do(0 ) 
ey 
Cancelling terms yields, 
+ 
: +, _ B do(0™) 
een 08) u dx 


Pp 


Rearranging, the Dankwerts' boundary condition at x = 0 is 


u 
erate = Fis — Pd) 32 8 
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From equation D.13, and noting that A, is a constant, we 


have 
oD = 0 
ax 
and 
do, (L") 
— 0 
ax 


Substituting these expressions and equation D.7 into equation 
Dec , 
Poag ee) (Lb) - =2.() 
P¢ Pp u ax an P¢ u 
Pp @ 
Upon cancelling terms, the second Dankwets' boundary condi- 


tion becomes, 


An important underlying consideration for using the 
Dankwerts' boundary conditions is that it simplifies the 
analysis since equation D.2 may be solved independently 


without haveing to consider entrance and exit regions. 


al. 





APPENDIX E 


DERIVATION OF THE FEM OPERATORS 


In the section on the finite element formulation, the 


following five differential operators were identified, 


1 

ero (E.1) 
1 

jG dG sean Ge) 
1 

ST at (E23) 
il 

i (G5 }<Gi>{¥}<G.>dn (E.4) 
0 
1 

J {G,}dn (E35) 
0 


where the G; are the global basis functions. These opera- 
tors are constructed on the element level by introducing the 
corresponding element basis functions, J; - The global and 


element basis functions are related by, 
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fee? 6h) 


where or and G5 are defined by 


(l - =) for — in element (e) 
(e) _ i 
2) a 


0 for &— not in element (e) 


for &— in element (e) 


o | 


3) 


0 for — not in element (e) 


and ae is the length of the (e)th element. 

The derivation of the local elemental matrices according 
to the Galerkin method for the global operations E.1 through 
E.5 proceeds as follows: 


For operator E.l, 


global local 


ab i (91 


J {Gj}<Gi>dn » uf <1 g5>dé] 
0 


Noting that, 
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Oe 





ee EEE 


bic ie 2 
Q cy (5) a ant 
e e 2 1 
, Ean’ Wiz 
ib ee Lien, 
a epee cee ab E 
e e 
mee Operator E.2, 
global local 
1 *e (91 
f (G4 }<Gj>dn » OU J <gi g3>dé] 
32 


Substituting in the local shape functions gives, 


E 
; (l - z. 
il 1 
j <= i (5—) >dé 
e e 


Carrying out the operations gives 
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i E ik 2 
a f° oo! oe 
e e e e 2 l 
J dé bp 5 
0 2 2 
- (5) (5) 
Pag) 22 
e 
Mor operator £.3, 
global local 
se 72 


aL 
SF ¥ cient <g, g5>dé] 


Substituting in for the local shape functions 


gE 
(1 - cam 
: Seats 
i SS ea) ace 
e e 


2g 
Cae 
e 


the elemental matrix becomes 


2 2 
Bats 2 + +5) | ar 4 
e e V2 e QR 
eS S e 
j dc ph > 
: E =e 2 ae 
as. Es 1-224 & 
de 2 Q a 
S (= 
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For operator E.4, 





global local 
1 te ai) y 
J co WSS Cat bY Oss <9195? <J195>dé} 
0 0 | 
35 bs 


Multiplying out the matrices, the expression becomes, 


| | f t 
Le Cao eee ade W495 9495*5-1 + 93,9995" ;) 
0 ae 
| f t 
(9597954%5_3 + 945999%3) $9 9990%3_3 + 5999T9";) 


Integrating each term separately, 


Q 
a gd = Se ne a = _i 
oc ai a Z Th a i 3 
0 0 e Le ae 
= es 2 
a0 Ze E 1 
| GC GCL eee 
0 1°12 0 2 ue u° 3 
des ae E z¢ 1 
fee Sr te = 
0 0 am es 
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Me se 


& & : re 1 
S gy9n9gdt = fo ~y > ty db = & 
0 0 2 g 
= = 
Ms + 
t = — — 
Q 
: 1 
f, 91929348 = & 
aa Me re 1 
| = ay Fs = = one 
J gigngnde = f x dé 3 
0 0 Q 
S 
‘ oa : 
G,g,9,d— = oe 5 
ae Ae 3 
= 
Substituting the values into the expression 
1 i 1 1 
eon eee 3° * 3) CE G-1)° 6 i’ 
1 i 1 1 
momo ree) 3 4-1 * 3°)? 


factoring out a - 1/3, the local matrix for the global 


operator becomes 
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eth ak 
, (ode ee) ( 5 ) 
1 
nee} <G.>1¥)<G.>d - = 
\ —_— all » 2 me 
LS aL 
mer eperator E.5, 
global local 
af “es cal 
f {G,jdn ul f dé] 
0 0 
32 


Substituting in the local shape function, and integrating, 


the expression becomes 


» (a2) 1 
S S ae 
j dp + 
0 
a i 
Q 
e 


This last operator is used for the excitation vector as 


described in the FEM formulation. 
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FIGURE 1 Semenov model of combustion. 








FIGURE 2 Idealized geometry of a fibrous graphite 
plate. 


99 











eeneeoe ve oe ¢ 8 & 2 o 
oon? 78878 © 8 @ 


e 


oo #0 @ 


oy 
e 
° 
e 
e 
* 


. 
oot 2 @ 8 8 @ ee 
eeteoeevweoee* ## © #8 
e eo? 
eet ee Pe ee 
eeeoeevov, eee sree ee ee 


‘ 
ofa oe? 0 0 0 Oo @ 





FIGURE 3 Separating a differential volume of porous medium 
into respective volumes of fiber and air. : 
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FIGURE 5 Energy balance on the internal flow. 
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FIGURE 6 Molecule-mass balance on the oxygen. 


———EEEEE 


| 100 








ONC > : Grr 


i 2 3 N-2 n-| n 
FIGURE 7 Linear shape functions for the Galerkin formulation 
of FEM. 
Element (1) Cay. a (n-2) (n-1) 
Local nodal point |! 2 2 > Si nn oe Daal 2 
Global nodal point 3 2 3 n-2 n-t 
6), Ma Q, We O, = om 
Vel VI, Wer. Yer, 
WED, YEP, WF, Ye, 


FIGURE 8 Numbering scheme used in system matrices. 
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FIGURE 9 Shows coupling in system of p. d. e. 
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